Data Visualization of WT and Nrf2 KO (-/-) BL6 PEITC or Cranberry Treated Mice 16S Microbiome Data Analysis, September 2019 Batch

Date: Thu Feb 27 13:12:41 2020
Scientist: Ran Yin
Sequencing (Waksman): Dibyendu Kumar
Statistics: Davit Sargsyan
Principal Investigator: Ah-Ng Kong

# Taxonomic Ranks:
# **K**ing **P**hillip **C**an n**O**t **F**ind **G**reen **S**ocks
# * Kingdom                
# * Phylum                    
# * Class                   
# * Order                   
# * Family     
# * Genus     
# * Species  

```r
# options(stringsAsFactors = FALSE,
#         scipen = 999)
# # Increase mmemory size to 64 Gb----
# invisible(utils::memory.limit(65536))
# str(knitr::opts_chunk$get())
# # NOTE: the below does not work!
# knitr::opts_chunk$set(echo = FALSE, 
#                       message = FALSE,
#                       warning = FALSE,
#                       error = FALSE)
# require(knitr)
# require(kableExtra)
require(phyloseq)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBwaHlsb3NlcVxuIn0= -->

Loading required package: phyloseq

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuIyByZXF1aXJlKHNoaW55KVxucmVxdWlyZShkYXRhLnRhYmxlKVxuYGBgXG5gYGAifQ== -->

```r
```r
# require(shiny)
require(data.table)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBkYXRhLnRhYmxlXG5kYXRhLnRhYmxlIDEuMTIuMiB1c2luZyAxOCB0aHJlYWRzIChzZWUgP2dldERUdGhyZWFkcykuICBMYXRlc3QgbmV3czogci1kYXRhdGFibGUuY29tXG4ifQ== -->

Loading required package: data.table data.table 1.12.2 using 18 threads (see ?getDTthreads). Latest news: r-datatable.com

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucmVxdWlyZShnZ3Bsb3QyKVxuYGBgXG5gYGAifQ== -->

```r
```r
require(ggplot2)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBnZ3Bsb3QyXG4ifQ== -->

Loading required package: ggplot2

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucmVxdWlyZShwbG90bHkpXG5gYGBcbmBgYCJ9 -->

```r
```r
require(plotly)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBwbG90bHlcblxuQXR0YWNoaW5nIHBhY2thZ2U6IMOi4oKsy5xwbG90bHnDouKCrOKEolxuXG5UaGUgZm9sbG93aW5nIG9iamVjdCBpcyBtYXNrZWQgZnJvbSDDouKCrMuccGFja2FnZTpnZ3Bsb3Qyw6LigqzihKI6XG5cbiAgICBsYXN0X3Bsb3RcblxuVGhlIGZvbGxvd2luZyBvYmplY3QgaXMgbWFza2VkIGZyb20gw6LigqzLnHBhY2thZ2U6c3RhdHPDouKCrOKEojpcblxuICAgIGZpbHRlclxuXG5UaGUgZm9sbG93aW5nIG9iamVjdCBpcyBtYXNrZWQgZnJvbSDDouKCrMuccGFja2FnZTpncmFwaGljc8Oi4oKs4oSiOlxuXG4gICAgbGF5b3V0XG4ifQ== -->

Loading required package: plotly

Attaching package: ‘plotly’

The following object is masked from ‘package:ggplot2’:

last_plot

The following object is masked from ‘package:stats’:

filter

The following object is masked from ‘package:graphics’:

layout

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucmVxdWlyZShEVClcbmBgYFxuYGBgIn0= -->

```r
```r
require(DT)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBEVFxuIn0= -->

Loading required package: DT

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxucmVxdWlyZShsbWVyVGVzdClcbmBgYFxuYGBgIn0= -->

```r
```r
require(lmerTest)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBsbWVyVGVzdFxuTG9hZGluZyByZXF1aXJlZCBwYWNrYWdlOiBsbWU0XG5Mb2FkaW5nIHJlcXVpcmVkIHBhY2thZ2U6IE1hdHJpeFxuXG5BdHRhY2hpbmcgcGFja2FnZTogw6LigqzLnGxtZXJUZXN0w6LigqzihKJcblxuVGhlIGZvbGxvd2luZyBvYmplY3QgaXMgbWFza2VkIGZyb20gw6LigqzLnHBhY2thZ2U6bG1lNMOi4oKs4oSiOlxuXG4gICAgbG1lclxuXG5UaGUgZm9sbG93aW5nIG9iamVjdCBpcyBtYXNrZWQgZnJvbSDDouKCrMuccGFja2FnZTpzdGF0c8Oi4oKs4oSiOlxuXG4gICAgc3RlcFxuIn0= -->

Loading required package: lmerTest Loading required package: lme4 Loading required package: Matrix

Attaching package: ‘lmerTest’

The following object is masked from ‘package:lme4’:

lmer

The following object is masked from ‘package:stats’:

step

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuc291cmNlKFxcc291cmNlL2Z1bmN0aW9uc19tYXkyMDE5LlJcXClcbiMgT24gV2luZG93cyBzZXQgbXVsdGl0aHJlYWQ9RkFMU0UtLS0tXG5tdCA8LSBUUlVFXG5gYGBcbmBgYCJ9 -->

```r
```r
source(\source/functions_may2019.R\)
# On Windows set multithread=FALSE----
mt <- TRUE

<!-- rnb-source-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


# Introduction
C57BL/6 wild-type (WT) and Nrf-2 double-knock-out (KO -/-) mice were given 2-week microbiome stabilization process using AIN93M diet and 8 more weeks to treat with either AIN93M or AIN93M 5% PEITC diet. Fecal samples were collected weekly, immediately frozen in liquid nitrogen and stored at -80^o^C. Serum, cecal, colon epithelial and whole colon tissues at week 10 were also collected for further analyses. Baseline, week 1 and 4 fecal samples were selected for 16s rRNA sequencing.  
  
This document examines results from the WT mice samples.  
  
We will attampt to answer the following questions:  
1. Did microbiome change over time?  
2. Was microbiome affected by diet?  
3. Was there a difference between the KO and WT?  
4. If there was a change in microbiome composition, what functional changes did it carry? What are the essential functions of the bacteria affected by the treatment and how can this be shown in vivo (metabolites, inflammation markers, etc.)?

# Data preprocessing
## Raw Data 
FastQ files were downloaded from [this Rutgers Box location](https://rutgers.app.box.com/folder/90143462291). A total of 144 files (2 per sample, pair-ended) and a pair of undetermined reads were downloaded. 

## Script
This script (***nrf2ubiome_dada2_sep2019_v1.Rmd***) was developed using [DADA2 Pipeline Tutorial (1.12)](https://benjjneb.github.io/dada2/tutorial.html) with tips and tricks from the [University of Maryland Shool of Medicine Institute for Genome Sciences (IGS)](http://www.igs.umaryland.edu/) [Microbiome Analysis Workshop (April 8-11, 2019)](http://www.igs.umaryland.edu/education/wkshp_metagenome.php). The output of the DADA2 script (***data_may2019/ps_sep2019.RData***) is explored in this document.

# Meta data: sample description

<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuIyBMb2FkIGRhdGEtLS0tXG4jIENvdW50c1xubG9hZChcImRhdGFfc2VwMjAxOS9wc19zZXAyMDE5LlJEYXRhXCIpXG5cbiMgVGF4b25vbXlcbmxvYWQoXCJkYXRhX3NlcDIwMTkvdGF4YS5SRGF0YVwiKVxudGF4YSA8LSBkYXRhLnRhYmxlKHNlcTE2cyA9IHJvd25hbWVzKHRheGEpLFxuICAgICAgICAgICAgICAgICAgIHRheGEpXG5gYGAifQ== -->

```r
# Load data----
# Counts
load("data_sep2019/ps_sep2019.RData")

# Taxonomy
load("data_sep2019/taxa.RData")
taxa <- data.table(seq16s = rownames(taxa),
                   taxa)

NOTE: correction to the meta-data! (11/15/2019)

correct_samples <- fread("data_sep2019/16s metadata Sep-2019.csv")
ps_sep2019@sam_data$DSS <- correct_samples$DSS

1 Samples

ps_sep2019@sam_data$Genotype_Week <- paste(ps_sep2019@sam_data$genotype,
                                           ps_sep2019@sam_data$time,
                                           sep = "_")
ps_sep2019@sam_data$ID <- factor(paste0(ps_sep2019@sam_data$mice_num,
                                        ps_sep2019@sam_data$cage))

ps_sep2019@sam_data$TREATMENT <- paste0(ps_sep2019@sam_data$DSS,
                                        ps_sep2019@sam_data$PEITC,
                                        ps_sep2019@sam_data$cranberry)
ps_sep2019@sam_data$TREATMENT <- factor(ps_sep2019@sam_data$TREATMENT,
                                        levels = c("000",
                                                   "100",
                                                   "110",
                                                   "101"),
                                        labels = c("Naive",
                                                   "DSS",
                                                   "DSS+PEITC",
                                                   "DSS+Cranberry"))

samples <- ps_sep2019@sam_data
datatable(samples,
          options = list(pageLength = nrow(samples)))

2 Prune data

The OTUs were mapped to Bacteria (96.07%), Eukaryota (2.95%) and Archea (0.03%) kingdoms, and 75 OTUs (0.95%) undefined.

The total of 7,867 unique sequences were found. Out of those, 7,558 were mapped to bacterial genomes.

dim(ps_sep2019@otu_table@.Data)

[1]   72 7867

# Remove OTU not mapped to Bacteria
ps0 <- subset_taxa(ps_sep2019, 
                   Kingdom == "Bacteria")
dim(ps0@otu_table@.Data)

[1]   72 7558

Out of the 7,558 OTUs 7,247 belonged to 12 Phyla. 311 of the OTUs (or 4.11% of bacterial OTUs) could not be mapped to a phylum.

t2 <- data.table(table(tax_table(ps0)[, "Phylum"],
                                  exclude = NULL))
t2$V1[is.na(t2$V1)] <- "Unknown"
setorder(t2, -N)
t2[, pct := N/sum(N)]
setorder(t2, -N)

colnames(t2) <- c("Phylum",
                  "Number of OTUs",
                  "Percent of OTUs")

datatable(t2,
          rownames = FALSE,
          caption = "Number of Bacterial OTUs by Phylum",
          class = "cell-border stripe",
          options = list(search = FALSE,
                         pageLength = nrow(t2))) %>%
  formatCurrency(columns = 2,
                 currency = "",
                 mark = ",",
                 digits = 0) %>%
  formatPercentage(columns = 3,
                   digits = 2)

3 OTU table (first 10 rows)

4 Total counts per sample (i.e. sequencing depth)

t1 <- colSums(otu[, 7:ncol(otu)])
t1 <- data.table(SAMPLE_NAME = names(t1),
                 Total = t1)

t2 <- data.table(SAMPLE_NAME = rownames(samples),
                 ID = samples$ID,
                 CAGE = samples$cage,
                 TREATMENT = samples$TREATMENT,
                 Genotype = samples$genotype,
                 WEEK = samples$time)

smpl <- merge(t1,
              t2,
              by = "SAMPLE_NAME")

p1 <- ggplot(smpl,
             aes(x = SAMPLE_NAME,
                 y = Total,
                 fill = TREATMENT,
                 colour = WEEK)) +
  facet_wrap(~ Genotype,
             scale = "free_x") +
  geom_bar(stat = "identity") +
  scale_x_discrete("") +
  scale_y_continuous("Number of Reads") +
  scale_fill_discrete("Treatment") +
  theme(axis.text.x = element_text(angle = 45,
                                   hjust = 1)) 
ggplotly(p1)

5 Richness (Alpha diversity)

Shannon index (aka Shannon enthrophy) is calculated as:
H’ = -sum(1 to R)p(i)ln(p(i)) When there is exactly 1 type of data (e.g. a single species in the sample), H’=0. The opposite scenario is when there are R>1 species present in the sample in the exact same amounts and H’=ln(R).

Shannon’s diversity index was calculated for each sample and ploted over time using the 7,764 from the 13 Phylum above.

shannon.ndx <- estimate_richness(ps0,
                                 measures = "Shannon")

shannon.ndx <- data.table(SAMPLE_NAME = rownames(shannon.ndx),
                          shannon.ndx)

smpl <- merge(smpl,
              shannon.ndx,
              by = "SAMPLE_NAME")

p1 <- ggplot(smpl,
             aes(x = Total,
                 y = Shannon,
                 fill = Genotype,
                 shape = WEEK)) +
  geom_point(size = 2) +
  scale_shape_manual(breaks = unique(smpl$WEEK),
                     values = 21:23)

tiff(filename = "tmp/shannon_vs_depth.tiff",
     height = 5,
     width = 6,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

ggplotly(p1)

Even though estimate_richness function does not adjust for the sequencing depth, there is no correlation between the index and the sample’s sequecing depth. Proceed with the comparison.

6 Shannon idex over time

p1 <- plot_richness(ps0,
                    x = "time", 
                    measures = "Shannon") +
  facet_wrap(~ genotype) +
  geom_line(aes(group = ID),
            color = "black") +
  geom_point(aes(fill = TREATMENT),
             shape = 21,
             size = 3,
             color = "black") +
  scale_x_discrete("") +
  theme(axis.text.x = element_text(angle = 30,
                                   hjust = 1,
                                   vjust = 1))

ggplotly(p = p1,
         tooltip = c("ID",
                     "value"))

p1 <- p1 + 
  scale_fill_discrete("") +
  theme(legend.position = "top")

tiff(filename = "tmp/shannon.tiff",
     height = 4,
     width = 5,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

The plot above suggests that the largest differences in alpha diversity (as measured by Shannon’s index) are in genotype.

Test if the richness changed between the baseline and Week 8.

smpl$TREATMENT <- factor(smpl$TREATMENT,
                         levels = c("DSS",
                                    "Naive",
                                    "DSS+PEITC",
                                    "DSS+Cranberry"))

tmp <- droplevels(smpl[WEEK != "week1"])

m1 <- lm(Shannon  ~ WEEK*(TREATMENT + Genotype),
         # offset = Total,
         data = tmp)
summary(m1)

Call:
lm(formula = Shannon ~ WEEK * (TREATMENT + Genotype), data = tmp)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.316186 -0.091027  0.007886  0.110704  0.293230 

Coefficients:
                                 Estimate Std. Error t value Pr(>|t|)
(Intercept)                       5.94987    0.07064  84.233  < 2e-16
WEEKweek8                         0.01158    0.09989   0.116   0.9084
TREATMENTNaive                    0.14581    0.08935   1.632   0.1109
TREATMENTDSS+PEITC               -0.03923    0.08935  -0.439   0.6631
TREATMENTDSS+Cranberry           -0.22582    0.08935  -2.527   0.0158
Genotypewidetype                 -0.54156    0.06318  -8.572 2.06e-10
WEEKweek8:TREATMENTNaive          0.01181    0.12636   0.093   0.9261
WEEKweek8:TREATMENTDSS+PEITC      0.01652    0.12636   0.131   0.8966
WEEKweek8:TREATMENTDSS+Cranberry  0.21535    0.12636   1.704   0.0965
WEEKweek8:Genotypewidetype        0.23085    0.08935   2.584   0.0137
                                    
(Intercept)                      ***
WEEKweek8                           
TREATMENTNaive                      
TREATMENTDSS+PEITC                  
TREATMENTDSS+Cranberry           *  
Genotypewidetype                 ***
WEEKweek8:TREATMENTNaive            
WEEKweek8:TREATMENTDSS+PEITC        
WEEKweek8:TREATMENTDSS+Cranberry .  
WEEKweek8:Genotypewidetype       *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1548 on 38 degrees of freedom
Multiple R-squared:  0.7845,    Adjusted R-squared:  0.7335 
F-statistic: 15.37 on 9 and 38 DF,  p-value: 3.671e-10

m2 <- lmer(Shannon  ~ WEEK*(TREATMENT + Genotype) + (1 | ID),
           # offset = Total,
           data = tmp)
summary(m2)

Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: Shannon ~ WEEK * (TREATMENT + Genotype) + (1 | ID)
   Data: tmp

REML criterion at convergence: -22.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.53721 -0.47495  0.06753  0.44489  1.53874 

Random effects:
 Groups   Name        Variance Std.Dev.
 ID       (Intercept) 0.01259  0.1122  
 Residual             0.01136  0.1066  
Number of obs: 48, groups:  ID, 24

Fixed effects:
                                 Estimate Std. Error       df t value
(Intercept)                       5.94987    0.07064 29.76837  84.233
WEEKweek8                         0.01158    0.06879 19.00000   0.168
TREATMENTNaive                    0.14581    0.08935 29.76837   1.632
TREATMENTDSS+PEITC               -0.03923    0.08935 29.76837  -0.439
TREATMENTDSS+Cranberry           -0.22582    0.08935 29.76837  -2.527
Genotypewidetype                 -0.54156    0.06318 29.76837  -8.572
WEEKweek8:TREATMENTNaive          0.01181    0.08701 19.00000   0.136
WEEKweek8:TREATMENTDSS+PEITC      0.01652    0.08701 19.00000   0.190
WEEKweek8:TREATMENTDSS+Cranberry  0.21535    0.08701 19.00000   2.475
WEEKweek8:Genotypewidetype        0.23085    0.06152 19.00000   3.752
                                 Pr(>|t|)    
(Intercept)                       < 2e-16 ***
WEEKweek8                         0.86814    
TREATMENTNaive                    0.11322    
TREATMENTDSS+PEITC                0.66379    
TREATMENTDSS+Cranberry            0.01704 *  
Genotypewidetype                 1.55e-09 ***
WEEKweek8:TREATMENTNaive          0.89350    
WEEKweek8:TREATMENTDSS+PEITC      0.85139    
WEEKweek8:TREATMENTDSS+Cranberry  0.02291 *  
WEEKweek8:Genotypewidetype        0.00135 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
                     (Intr) WEEKw8 TREATMENTN TREATMENTDSS+P
WEEKweek8            -0.487                                 
TREATMENTNv          -0.632  0.308                          
TREATMENTDSS+P       -0.632  0.308  0.500                   
TREATMENTDSS+C       -0.632  0.308  0.500      0.500        
Gentypwdtyp          -0.447  0.218  0.000      0.000        
WEEK8:TREATMENTN      0.308 -0.632 -0.487     -0.243        
WEEK8:TREATMENTDSS+P  0.308 -0.632 -0.243     -0.487        
WEEK8:TREATMENTDSS+C  0.308 -0.632 -0.243     -0.243        
WEEKwk8:Gnt           0.218 -0.447  0.000      0.000        
                     TREATMENTDSS+C Gntypw WEEK8:TREATMENTN
WEEKweek8                                                  
TREATMENTNv                                                
TREATMENTDSS+P                                             
TREATMENTDSS+C                                             
Gentypwdtyp           0.000                                
WEEK8:TREATMENTN     -0.243          0.000                 
WEEK8:TREATMENTDSS+P -0.243          0.000  0.500          
WEEK8:TREATMENTDSS+C -0.487          0.000  0.500          
WEEKwk8:Gnt           0.000         -0.487  0.000          
                     WEEK8:TREATMENTDSS+P WEEK8:TREATMENTDSS+C
WEEKweek8                                                     
TREATMENTNv                                                   
TREATMENTDSS+P                                                
TREATMENTDSS+C                                                
Gentypwdtyp                                                   
WEEK8:TREATMENTN                                              
WEEK8:TREATMENTDSS+P                                          
WEEK8:TREATMENTDSS+C  0.500                                   
WEEKwk8:Gnt           0.000                0.000              

7 Calculate change in Shannon index from baseline

dd <- smpl
dd[, delta := Shannon - Shannon[WEEK == "baseline"],
   by = ID]
dd$diff <- paste(dd$WEEK,
                 "-baseline",
                 sep = "")

dd <- dd[WEEK != "baseline",]

p1 <- ggplot(dd,
             aes(x = TREATMENT,
                 y = delta,
                 fill = Genotype)) +
  facet_wrap(~ diff) +
  geom_hline(yintercept = 0,
             linetype = "dashed") +
  geom_point(position = position_dodge(0.3),
             shape = 21,
             size = 3) +
  scale_y_continuous("Shannon Index Percent Change from Baseline") +
  theme(axis.text.x = element_text(angle = 45,
                                   hjust = 1))
print(p1)

dd$TREATMENT <- factor(dd$TREATMENT,
                        levels = c("DSS",
                                   "Naive",
                                   "DSS+PEITC",
                                   "DSS+Cranberry"))
dd$Genotype <- factor(dd$Genotype,
                       levels = c("widetype",
                                  "nrf2KO"))

m1 <- lm(delta ~ TREATMENT*Genotype,
         data = dd)
summary(m1)

Call:
lm(formula = delta ~ TREATMENT * Genotype, data = dd)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.40513 -0.09560 -0.02012  0.09568  0.35517 

Coefficients:
                                      Estimate Std. Error t value
(Intercept)                            0.25142    0.07286   3.451
TREATMENTNaive                        -0.04426    0.10303  -0.430
TREATMENTDSS+PEITC                    -0.15777    0.10303  -1.531
TREATMENTDSS+Cranberry                 0.04463    0.10303   0.433
Genotypenrf2KO                        -0.18412    0.10303  -1.787
TREATMENTNaive:Genotypenrf2KO         -0.02851    0.14571  -0.196
TREATMENTDSS+PEITC:Genotypenrf2KO      0.24747    0.14571   1.698
TREATMENTDSS+Cranberry:Genotypenrf2KO  0.07927    0.14571   0.544
                                      Pr(>|t|)   
(Intercept)                            0.00133 **
TREATMENTNaive                         0.66985   
TREATMENTDSS+PEITC                     0.13358   
TREATMENTDSS+Cranberry                 0.66720   
Genotypenrf2KO                         0.08153 . 
TREATMENTNaive:Genotypenrf2KO          0.84586   
TREATMENTDSS+PEITC:Genotypenrf2KO      0.09721 . 
TREATMENTDSS+Cranberry:Genotypenrf2KO  0.58946   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1785 on 40 degrees of freedom
Multiple R-squared:  0.249, Adjusted R-squared:  0.1176 
F-statistic: 1.894 on 7 and 40 DF,  p-value: 0.09608

# No significant interactions, proceed with 2-way analysis
m2 <- lm(delta ~ TREATMENT + Genotype,
         data = dd)
summary(m2)

Call:
lm(formula = delta ~ TREATMENT + Genotype, data = dd)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.49158 -0.09742 -0.01290  0.11101  0.35281 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)             0.21414    0.05849   3.661 0.000683 ***
TREATMENTNaive         -0.05851    0.07399  -0.791 0.433377    
TREATMENTDSS+PEITC     -0.03404    0.07399  -0.460 0.647813    
TREATMENTDSS+Cranberry  0.08427    0.07399   1.139 0.261014    
Genotypenrf2KO         -0.10956    0.05232  -2.094 0.042176 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1812 on 43 degrees of freedom
Multiple R-squared:  0.1674,    Adjusted R-squared:  0.08999 
F-statistic: 2.162 on 4 and 43 DF,  p-value: 0.08947

At Week 8 there was significantly smaller increase of alpha diversity from baseline in Nrf2 KO compared to WT, and in DSS+Cranberry compared to DSS only.

8 Load aminoacids

aa <- fread("data_sep2019/sep2019_aminoacids.csv")

aa <- aa[!is.na(ID), ]
aa$ID <- paste0(aa$ID,
                aa$CAGE)

smpl1 <- unique(smpl[, c("ID",
                            "TREATMENT",
                            "Genotype")])
smpl1$ID <- as.character(smpl1$ID)
aa <- merge(smpl1,
            aa,
            by = "ID")
aa[, trt_week := paste(TREATMENT,
                       WEEK,
                       sep = "_")]
aa$trt_week <- factor(aa$trt_week,
                      levels = c("Naive_week2",
                                 "Naive_week6",
                                 "DSS_week2",
                                 "DSS_week6",
                                 "DSS+Cranberry_week2",
                                 "DSS+Cranberry_week6" ,
                                 "DSS+PEITC_week2",
                                 "DSS+PEITC_week6"))

for (i in 8:(ncol(aa) - 1)) {
  tmp <- aa[, c(1, 3, 7, ncol(aa), i), with = FALSE]
  colnames(tmp)[5] <- "Y"
  p1 <- ggplot(tmp,
               aes(x = trt_week,
                   y = Y,
                   fill = Genotype,
                   group = ID)) +
    geom_line(position = position_dodge(0.3)) +
    geom_point(shape = 21,
               size = 3,
               position = position_dodge(0.3)) +
    scale_x_discrete("") +
    scale_y_continuous(colnames(aa)[i]) +
    theme(axis.text.x = element_text(angle = 45,
                                     hjust = 1))
  # tiff(filename = paste0("tmp/",
  #                        colnames(aa)[i],
  #                        ".tiff"),
  #      height = 4,
  #      width = 5,
  #      units = "in",
  #      res = 600,
  #      compression = "lzw+p")
  # print(p1)
  # graphics.off()
  
  print(p1)
}

9 Aminoacid data PCA

dt_pca <- aa[, Alanine:glutamine]

m1 <- prcomp(dt_pca)
summary(m1)

Importance of components:
                          PC1    PC2     PC3     PC4     PC5     PC6
Standard deviation     4.8108 2.5833 1.78182 1.13495 0.93367 0.84949
Proportion of Variance 0.6276 0.1810 0.08609 0.03493 0.02364 0.01957
Cumulative Proportion  0.6276 0.8086 0.89466 0.92959 0.95323 0.97280
                           PC7     PC8     PC9    PC10    PC11    PC12
Standard deviation     0.62035 0.49052 0.35960 0.26786 0.22940 0.20322
Proportion of Variance 0.01044 0.00652 0.00351 0.00195 0.00143 0.00112
Cumulative Proportion  0.98323 0.98976 0.99326 0.99521 0.99663 0.99775
                         PC13    PC14    PC15    PC16    PC17    PC18
Standard deviation     0.1486 0.13250 0.11932 0.11154 0.08896 0.07274
Proportion of Variance 0.0006 0.00048 0.00039 0.00034 0.00021 0.00014
Cumulative Proportion  0.9983 0.99883 0.99922 0.99955 0.99977 0.99991
                          PC19    PC20
Standard deviation     0.05110 0.02565
Proportion of Variance 0.00007 0.00002
Cumulative Proportion  0.99998 1.00000

# Select PC-s to pliot (PC1 & PC2)
choices <- 1:2
# Scores, i.e. points (df.u)
dt.scr <- data.table(m1$x[, choices])
# Add grouping variable
dt.scr$grp <- aa$trt_week
dt.scr$TREATMENT <- aa$TREATMENT
dt.scr$WEEK <- aa$WEEK
dt.scr

# Loadings, i.e. arrows (df.v)
dt.rot <- as.data.frame(m1$rotation[, choices])
dt.rot$feat <- rownames(dt.rot)
dt.rot <- data.table(dt.rot)
dt.rot

dt.load <- melt.data.table(dt.rot,
                           id.vars = "feat",
                           measure.vars = 1:2,
                           variable.name = "pc",
                           value.name = "loading")
dt.load$feat <- factor(dt.load$feat,
                       levels = unique(dt.load$feat))
# Plot loadings
p0 <- ggplot(data = dt.load,
             aes(x = feat,
                 y = loading)) +
  facet_wrap(~ pc,
             nrow = 2) +
  geom_bar(stat = "identity") +
  ggtitle("PC Loadings") +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))
tiff(filename = "tmp/pc.1.2_loadings.tiff",
     height = 5,
     width = 8,
     units = 'in',
     res = 300,
     compression = "lzw+p")
print(p0)
graphics.off()

print(p0)

# Axis labels
u.axis.labs <- paste(colnames(dt.rot)[1:2], 
                     sprintf('(%0.1f%% explained var.)', 
                             100*m1$sdev[choices]^2/sum(m1$sdev^2)))
u.axis.labs

[1] "PC1 (62.8% explained var.)" "PC2 (18.1% explained var.)"

# Based on Figure p0, keep only a few variables with high loadings in PC1 and PC2----
# var.keep.ndx <- which(dt.rot$feat %in% c(...))
# Or select all
# var.keep.ndx <- 3:ncol(dt1)
# Use dt.rot[var.keep.ndx,] and dt.rot$feat[var.keep.ndx]
p1 <- ggplot(data = dt.rot,
             aes(x = PC1,
                 y = PC2)) +
  coord_equal() +
  geom_point(data = dt.scr,
             aes(fill = grp),
             shape = 21,
             size = 2,
             alpha = 0.5) +
  geom_segment(aes(x = 0,
                   y = 0,
                   xend = 10*PC1,
                   yend = 10*PC2),
               arrow = arrow(length = unit(1/2, 'picas')),
               size = 1.2, 
               color = "black") +
  geom_text(aes(x = 11*PC1,
                y = 11*PC2,
                label = dt.rot$feat),
            # size = 5,
            hjust = 0.5) +
  scale_x_continuous(u.axis.labs[1]) +
  scale_y_continuous(u.axis.labs[2]) +
  scale_fill_discrete(name = "Group") +
  ggtitle("Biplot of Aminoacids") +
  theme(plot.title = element_text(hjust = 0.5,
                                  size = 20))
tiff(filename = "tmp/aminoacids_biplot.tiff",
     height = 10,
     width = 10,
     units = 'in',
     res = 300,
     compression = "lzw+p")
print(p1)
graphics.off()

ggplotly(p1)

p2 <- ggplot(data = dt.rot,
             aes(x = PC1,
                 y = PC2)) +
  coord_equal() +
  geom_point(data = dt.scr,
             aes(fill = WEEK),
             shape = 21,
             size = 2,
             alpha = 0.5) +
  geom_segment(aes(x = 0,
                   y = 0,
                   xend = 10*PC1,
                   yend = 10*PC2),
               arrow = arrow(length = unit(1/2, 'picas')),
               size = 1.2, 
               color = "black") +
  geom_text(aes(x = 11*PC1,
                y = 11*PC2,
                label = dt.rot$feat),
            # size = 5,
            hjust = 0.5) +
  scale_x_continuous(u.axis.labs[1]) +
  scale_y_continuous(u.axis.labs[2]) +
  scale_fill_discrete(name = "Week") +
  ggtitle("Biplot of Aminoacids") +
  theme(plot.title = element_text(hjust = 0.5,
                                  size = 20))
tiff(filename = "tmp/aminoacids_by_week_biplot.tiff",
     height = 10,
     width = 10,
     units = 'in',
     res = 300,
     compression = "lzw+p")
print(p2)
graphics.off()

ggplotly(p2)

p2 <- ggplot(data = dt.rot,
             aes(x = PC1,
                 y = PC2)) +
  coord_equal() +
  geom_point(data = dt.scr,
             aes(fill = TREATMENT),
             shape = 21,
             size = 2,
             alpha = 0.5) +
  geom_segment(aes(x = 0,
                   y = 0,
                   xend = 10*PC1,
                   yend = 10*PC2),
               arrow = arrow(length = unit(1/2, 'picas')),
               size = 1.2, 
               color = "black") +
  geom_text(aes(x = 11*PC1,
                y = 11*PC2,
                label = dt.rot$feat),
            # size = 5,
            hjust = 0.5) +
  scale_x_continuous(u.axis.labs[1]) +
  scale_y_continuous(u.axis.labs[2]) +
  scale_fill_discrete(name = "Treatment") +
  ggtitle("Biplot of Aminoacids") +
  theme(plot.title = element_text(hjust = 0.5,
                                  size = 20))
tiff(filename = "tmp/aminoacids_by_trt_biplot.tiff",
     height = 10,
     width = 10,
     units = 'in',
     res = 300,
     compression = "lzw+p")
print(p2)
graphics.off()

ggplotly(p2)

10 Remove unmapped OTUs

The 311 unmapped OTUs were removed from further analysis (with 7,247 OTUs left).

ps1 <- subset_taxa(ps0, 
                   !is.na(Phylum))
dim(ps1@otu_table@.Data)

[1]   72 7247

11 Counts at Phylum level

12 Relative abundance (%) at Phylum level

Remove phyla with relative abundance of >= 1% in less than 10% of samples.

t1 <- data.table(Phylum = ra_p$Phylum,
                 `Number of Samples` = rowSums(ra_p[, 2:ncol(ra_p)] >= 0.01))
t1$`Percent Samples` <-  t1$`Number of Samples`/72

setorder(t1, -`Number of Samples`)
datatable(t1,
          rownames = FALSE,
          caption = "Taxonomic  count table",
          class = "cell-border stripe",
          options = list(search = FALSE,
                         pageLength = nrow(t1))) %>%
  formatPercentage(columns = 3,
                   digits = 1)

We will remove Chlamydiae from this analysis.

XX OTUs, down from YY OTUs in the previous table.

13 Relative Abundance in Samples at Different Taxonomic Ranks

13.1 1. Class

mu$Trt_Genotype <- factor(paste(mu$Treatment,
                          mu$Genotype,
                          sep = "_"))

p0 <- ggplot(mu,
             aes(x = Week,
                 y = x,
                 group = Trt_Genotype)) +
  facet_wrap(~ Class,
             scale = "free_y") +
  geom_line(position = position_dodge(0.3)) +
  geom_point(aes(fill = Trt_Genotype),
             shape = 21,
             size = 2,
             alpha = 0.5,
             position = position_dodge(0.3)) +
  scale_x_discrete("") +
  scale_y_continuous("Relative Abundance (%)") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))

tiff(filename = "tmp/wt_class_over_time.tiff",
     height = 5,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p0)
graphics.off()

print(p0)

p1 <- ggplot(mu,
             aes(x = x,
                 y = Class,
                 color = Trt_Genotype,
                 shape = Week)) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)") +
  theme(legend.position = "top")

tiff(filename = "tmp/wt_class_ra.tiff",
     height = 4,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

ggplotly(p1)

13.2 2. Order

mu$Trt_Genotype <- factor(paste(mu$Treatment,
                                mu$Genotype,
                                sep = "_"))

p0 <- ggplot(mu,
             aes(x = Week,
                 y = x,
                 group = Trt_Genotype)) +
  facet_wrap(~ Order,
             scale = "free_y") +
  geom_line(position = position_dodge(0.3)) +
  geom_point(aes(fill = Trt_Genotype),
             shape = 21,
             size = 2,
             alpha = 0.5,
             position = position_dodge(0.3)) +
  scale_x_discrete("") +
  scale_y_continuous("Relative Abundance (%)") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))

tiff(filename = "tmp/wt_Order_over_time.tiff",
     height = 5,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p0)
graphics.off()

print(p0)

p1 <- ggplot(mu,
             aes(x = x,
                 y = Order,
                 color = Trt_Genotype,
                 shape = Week)) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)") +
  theme(legend.position = "top")

tiff(filename = "tmp/wt_Order_ra.tiff",
     height = 4,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

ggplotly(p1)

13.3 3. Family

NOTE: only the first 24 families had large enough counts - ploting only them.

mu$Trt_Genotype <- factor(paste(mu$Treatment,
                                mu$Genotype,
                                sep = "_"))
mu1 <- droplevels(mu[Family %in% levels(mu$Family)[nlevels(mu$Family):(nlevels(mu$Family) - 24)], ])

p0 <- ggplot(mu1,
             aes(x = Week,
                 y = x,
                 group = Trt_Genotype)) +
  facet_wrap(~ Family,
             scale = "free_y") +
  geom_line(position = position_dodge(0.3)) +
  geom_point(aes(fill = Trt_Genotype),
             shape = 21,
             size = 2,
             alpha = 0.5,
             position = position_dodge(0.3)) +
  scale_x_discrete("") +
  scale_y_continuous("Relative Abundance (%)") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))

tiff(filename = "tmp/wt_Family_over_time.tiff",
     height = 7,
     width = 9,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p0)
graphics.off()

print(p0)

p1 <- ggplot(mu1,
             aes(x = x,
                 y = Family,
                 color = Trt_Genotype,
                 shape = Week)) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)") +
  theme(legend.position = "top")

tiff(filename = "tmp/wt_Family_ra.tiff",
     height = 4,
     width = 7,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

ggplotly(p1)

13.4 4. Genus

mu$Trt_Genotype <- factor(paste(mu$Treatment,
                                mu$Genotype,
                                sep = "_"))
mu1 <- droplevels(mu[Genus %in% levels(mu$Genus)[nlevels(mu$Genus):(nlevels(mu$Genus) - 35)], ])

p0 <- ggplot(mu1,
             aes(x = Week,
                 y = x,
                 group = Trt_Genotype)) +
  facet_wrap(~ Genus,
             scale = "free_y") +
  geom_line(position = position_dodge(0.3)) +
  geom_point(aes(fill = Trt_Genotype),
             shape = 21,
             size = 2,
             alpha = 0.5,
             position = position_dodge(0.3)) +
  scale_x_discrete("") +
  scale_y_continuous("Relative Abundance (%)") +
  theme(legend.position = "top",
        axis.text.x = element_text(angle = 45,
                                   hjust = 1))

tiff(filename = "tmp/wt_Genus_over_time.tiff",
     height = 9,
     width = 12,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p0)
graphics.off()

print(p0)

p1 <- ggplot(mu1,
             aes(x = x,
                 y = Genus,
                 color = Trt_Genotype,
                 shape = Week)) +
  geom_point(size = 3,
             alpha = 0.5) +
  geom_vline(xintercept = 1,
             linetype = "dashed") +
  scale_x_continuous("Relative Abundance (%)") +
  theme(legend.position = "top")

tiff(filename = "tmp/wt_Genus_ra.tiff",
     height = 9,
     width = 9,
     units = "in",
     res = 600,
     compression = "lzw+p")
print(p1)
graphics.off()

ggplotly(p1)

14 Session Information

sessionInfo()